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Abstract 

We investigate the collective phase dynamics in conventional long Josephson 
junction (LJJ) stacks and in layered superconductors, exhibiting intrinsic LJJ 
behaviors. Using a theoretical model which accounts for both the magnetic 
induction effect and the breakdown of local charge neutrality (i.e., charging 
effect), we show that the collective motion of Josephson vortices, including 
the dispersion of Josephson plasma mode and the Swihart-type velocity, in an 
intrinsic LJJ stack such as b^S^CaC^Os+y (BSCCO) is significantly mod- 
ified from those in a conventional LJJ stack. In BSCCO, the strength of the 
charging effect a is small (i.e., a ~ 0.1 — 0.4), but it leads to notable changes in 
collective phase dynamics, including changes to the stability condition. Also, 
we show that splitting of the supercurrent branch in the resistive state is 
due to collective motion of Josephson vortices. The width of spread of these 
sub-branches in the linear current-voltage regime depends on a, suggesting 
another way to measure the charging effect in BSCCO. 
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I. INTRODUCTION 



Dynamics of magnetic vortices in a stack of long Josephson junctions (LJJ) in a magnetic 
field applied parallel to the junction layers have attracted much attention due to their 
intriguing applied and fundamental interests. The motion of Josephson vortices in a single 
junction system has been exploited in various devices.! Collective motion of these vortices in 
a LJJ stack, in which layers of superconductor (S) and insulator (I) are arranged vertically 
as in Fig. 1, can be exploited in high frequency devices! such as tunable submillimeter-wave 
oscillators and detectors. Here collective motion, including both in-phase and out-of-phase 
modes shown in Fig. 2, arises from mutual phase-locking of Josephson junctions caused by 
(magnetic) inductive coupling between screening currents flowing around adjacent Josephson 
vortices as they move under a bias current. The phase- locking establishes phase coherence 
across the Josephson junctions. A LJJ stack exhibiting this phase coherence leads to high 
power output and bandwidth, and it can serve as a model system for scientific studies.!'! 

The motion of Josephson vortices in LJJ stacks yields interesting phenomena: (i) 
Josephson plasma resonance!'! (JPR) and (ii) supercurrent sub-branching.011! The exper- 
iments on both conventional LJJ stacks! (e.g., Nb-Al/AlO^-Nb multilayers) and layered 
superconductors!®! (e.g., Bi 2 Sr 2 CaCu 2 8+ j / (BSCCO)) behaving as intrinsic LJJ stacks^ 
indicate that JPR can be tuned0 by magnetic field B and can occur over a broad range of 
frequencies, from microwave to submillimeter-wave. Also, the supercurrent branch in the 
current-voltage (I-V) data splits into multiple sub-branches when a bias current exceeds 
some critical value.! To explain the datajl!'!!'! two theoretical models have been proposed: 
one is based on the inductive coupling (i.e., magnetic induction model), and the other is 
based on the coupling due to charge variation in the S layers (i.e., charging effect model). 00 

The magnetic induction model assumes that the S layer thickness ds is much larger 
than the Debye (charge screening) length rp (i.e., ds ^> td), as in conventional LJJ. In 
this case, charge variations (or electric field) at each S layers are screen out, yielding local 
charge neutrality. Consequently, the electric field does not lead to the longitudinal coupling 



2 



between the S layers. In this model, an applied magnetic field induces supercurrents along 
the S layers and results in the inductive interactionEl'ElIll between adjacent S layers. The 
induction coupling strength is inversely proportional to the common S layer thickness. This 
model has been used to explain the experimental data for BSCCO. However, the underlying 
assumption is not justified in BSCCO since ds ~ 3A and r D ~ 2 — 3A.0'Ei 

On the other hand, the charging effect model accounts for the nonequilibrium effect in 
atomic scale thick superconducting layers. When the S layers are so thin to be comparable 
to the Debye length (i.e., ds ~ ro), as in BSCCO, the breakdown of local charge neutrality 
yields the charging effect. El The particle-hole imbalance^ may also occur since each super- 
conducting layers cannot completely screen out the charge variation. Hence the presence of 
charge variations yields the interaction between the contiguous superconducting layers and 
leads to the coupling between the S layers. Recently the charging effect model, neglecting 
the magnetic induction effect, has been used to interpret the data for BSCCO.ll! 

Earlier studies,00 including numerical simulations^^ of a finite L JJ system, show that 
these two models can explain the data qualitatively, but considerable inconsistencies between 
the experimental and the theoretical results have been found. For example, transverse and 
longitudinal JPR are predicted by the magnetic induction model and the charging effect 
model, respectively, but the data indicate that both types of resonance occur.0 Recent 
experimentsi'i on Hgl 2 -intercalated BSCCO and BSCCO single crystals indicate that the 
supercurrent branch in the I-V data splits into multiple sub-branches in the resistive state 
when B ~ H (i.e., low vortex density regime). An estimated valueilof H for Nb- Al/AlO^- 
Nb multilayers and BSCCO is roughly 0.001T and 0.2T, respectively. In the dense vortex 
regime (i.e., B 3> H a ), the I-V data exhibit characteristic kinks, and these kinks closely 
resemble the prediction made by Machida et al., using the magnetic induction model.il 
However a closer examinationi of the data reveals some inconsistencies.i These inconsisten- 
cies suggest that a better theoretical model is needed to describe the LJJ stacks. 

In this paper, we investigate collective phase dynamics in conventional LJJ stacks and 
layered superconductors at low magnetic fields (i.e., B ~ H Q in which Josephson vortices 



are in every I layers as in Fig. 2) and at low temperatures (i.e., below the Abrikosov vortex 
lattice melting temperature), using a theoretical model accounting for both the induction 
effect and the charging effect. These two effects are equally importantB^I in BSCCO since 
fD ~ ds-, but the charging effect is neglected in many studies because its strength a is small@ 
(e.g., a ~ 0.1 — 0.4 in BSCCO). We show how the collective motion of Josephson vortices is 
modified by a weak charging effect. We outline two main results. First, the Josephson plasma 
dispersion relation, the Swihart velocity, and the stability condition for collective motion 
in BSCCO are considerably modified from those in Nb-Al/AlO x -Nb multilayers. Second, 
the splitting of the supercurrent branch in the resistive state is due to collective motion of 
Josephson vortices, and the width of spread of these sub-branches in the linear I-V regime 
depends on a. These results are consistent with the experimental data described above. 

The remainder of the paper is organized as follows. In Sec. II, a theoretical model, 
which accounts for both the magnetic induction effect and the charging effect, is derived 
by extending previous models. In Sec. Ill, the Josephson plasma dispersion relation and 
the Swihart velocity for the collective modes are computed from our model derived in Sec. 
II. In Sec. IV, we determine the stability condition for the mutually phase-locked modes, 
performing the linear stability analysis. In Sec. V, we show that the splitting of the super- 
current branch in the resistive state is due to the collective motion of Josephson vortices. 
Finally, in Sec. VI, we summarize our results and conclude. 



In this section, we derive a theoretical model, extending previous approaches. A brief 
discussion of this model was published.il Here we consider a system with a large number of 
LJJ (i.e., N 3> 1) neglecting the boundary effect and present new results obtained from this 
model in later sections. To account for both the magnetic induction effect and the charging 
effect, we start with the gauge-invariant phase difference between the S layers I and i — 1, 



II. THEORETICAL MODEL 
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where 9t is the phase of the superconducting order parameter, <p = hc/2e is the flux quan- 
tum, and A is the vector potential in the I layers. In this paper, we employ the Cartesian 
coordinates and assume that the S and I layers are stacked along z-direction and the mag- 
netic field is applied along the ^/-direction, as in Fig 1. For simplicity, the thicknesses of the 
S (ds) and I (dj) layers are taken to be uniform. 

The magnetic induction effect due to the applied magnetic field (along the ^-direction) 
yields a spatial variation of the phase difference (along the ^-direction). An equation de- 
scribing the magnetic inductive couplings! between the S layers 

^ 1 = s(Bt + ij + JB£_x^_ 2 ) + d'Bij-x (2) 

is easily obtained by taking a spatial derivative of Eq. ([!]) and by using the expression for 
the supercurrent density 

V^-^AJ . (3) 



8tt 2 A 2 

Here d! = di + 2X coth(ds/A) and s = — X[sinh(ds / X)]~ l are expressed in terms of the London 
penetration depth A. The magnetic field -B^_i in the I layer between two S layers i and £ — 1 
is parallel to the layers. Note that B^_\ differs from B since the magnetic field generated by 
the supercurrent in the S layers modifies the field in the I layer. Using Maxwell's equation, 
we express the spatial derivative of the magnetic field as 

dB^ x 4tt 

q = —{Jc®n<Pt4-i-JB + Jtj-i) (4) 

where J c is the Josephson critical current density, and Jb is a bias current density. Note 
that the magnetic field entering the I layers yields0 a triangular Josephson vortex lattice 
(JVL) when the bias current is either absent or small. The current density^ J T , 

t t ^o_^_ d<Pt,e-i I dVtj-i f v 

J ^- 1 ~ 2nV dt + AttV dt ' [b) 

includes the quasiparticle and the displacement current contribution. Here T> = d'+2s = dj+ 
2A tanh(d,s/2A) is the effective thickness of the block layer, a is the quasiparticle conductivity, 



e is the dielectric constant of I layers, and V^_i is the voltage between the S layers i and 
£ — 1. Using Eq. (HI), we rewrite Eq. (0) as 



fr <9 2 c 



s(sin + sin ^_i / _ 2 ) + ci' sin . 



87T 2 <9x 2 

J' 7T 1 „ { jT 1 jT 



-VJ 



B 



+ ^ Jtji-i + + Ji-1,1-2) ■ (6) 

Note that S = s/d! measures the induction coupling strength, and S = —0.5 in the strong 
coupling limit. The phase difference equation (Eq. (4) in Ref. 27) derived by Bulaevskii and 
Clem within the framework of Lawrence-Doniach mo de@ can be obtained from Eq. @ when 
the time-dependent terms are neglected (i.e., Jj t _ x = 0) and the relations (87r 2 /0 o ) J c d' = 
(2/Aj) + (1/A 2 ) and (Sn 2 /(p )J c s = — 1/A 2 are used. Here A c is the magnetic penetration 
depth in the direction perpendicular to the S layers. 

The presence of a nonequilibrium state leads to the interaction between the S layers. 
When the S layer thickness is comparable to the Debye screening length (i.e., Td ~ ds), the 
S layers are in a nonequilibrium state because the charge variations in these layers are not 
completely screened. This incomplete charge screening enhances the temporal variation of 
the phase difference. One can include this effect in the phase dynamics, mo difyingHIll the 
usual AC Josephson relation, which is a time derivative of Eq. (P, to 
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Vtj-i + <&z- $ t -i (7) 



2tt dt 

as a way to account for a nonzero gauge-invariant potential = <f)g + (H/2e)(d9e/dt) gener- 
ated inside the S layers. Here <p£ is the electrostatic potential. The modified AC Josephson 
relation of Eq. (0) can be rewritten as 

<^C^£_l = _ a _ 2 + y J _y t+ ^ ^ (g) 

Z7T Ot 

using the charge density pi = — ($£ — \Ev)/47rrf> and the Maxwell's equation eV ■ E = 47rp. 
Here a = er 2 D jV>ds measures the strength of the charging effect, and measures the 
particle-hole imbalance in the S layer. For simplicity, we consider only the charging effect 
by setting = = 0, as it has been done in earlier studies.0 Note that the usual AC 
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Josephson relation is obtained from Eq. (|j) when a = 0. This indicates that the charging 
effect (i.e., a^O) enhances the coupling between neighboring junctions. Using Eq. (^), we 
relate the time derivative of the phase difference to the current densities and obtain 



Jr. 



1 <9 2 <^,£-i (3 dtpij-i a/3 d 

H ?n 777 Ve+ii - *<pe,t-i + 



v uj 2 dt 2 ujp dt uj p dt 



. rj~r rj~ \ 

Je,e-i ~ a [Je+i,t ~ 2Jg,e-i + Je-14-2) ■ (9) 

Here u p = c/(-y/eA c ) is the plasma frequency, (3 = (47r/c)(crA c /y / e) = and (3 C is the 

McCumber parameter. Note that the spatial variation of </^,i-i is neglected in the charging 
effect modelEl of Eq. (||). The terms of the order 0(a(3) can be safely neglected since a and 
(3 are small in the layered superconductors (i.e., a(3 <C 1). For example, the experimental 
value for a and (3 in BSCCO are roug hly 0.1-0.4 and 0.2, respectively.ffl Neglecting these 
small terms, we rewrite Eq. (|]) as 

(1 d ''(Pl i—\ (3 l rp rp rp rp 

TJI Qf2 ^ 7J J ~ Jt,i-i ~~ a \ Ji+i,i ~ ^Ji,e-i + J £-1,1-1) ■ (10) 

As we shall see in Sec. Ill, the charging effect terms in Eq. (|T0|) yield purely longitudinal 
Josephson plasma excitations Jil 

A theoretical model, including both the magnetic induction effect and the charging effect, 
can be obtained easily by noting that the magnetic induction model of Eq. @ and the 
charging effect model of Eq. fllCf ) are coupled to each other via the current density J T . 
Combining Eqs. (J6|) and fllQl) , we obtain the coupled sine-Gordon equations, 

1 ( 1 d 2 A e (3 dA e \ a 1 d 2 E e 



dx 2 \ 2 C V \lo 2 dt 2 uj p dt J X 2 c Vu 2 dt 2 
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X 2 C V 



d'sincpej-i + s(sin.cpe + i t i + sin^_i^_ 2 ) - V^j- 

'J c. 



(11) 



where A £ = d'ipi j e-i + s(<pt +lj £ + <^_i,*- 2 ) and E e = s(ip^ 2 ,£-3 + <f£+2,£+i) + (d'- 2s) (<^_i^_ 2 + 
<Pt+i,t) + 2(s — d')(f£ t £_i. The third term on the left hand side of Eq. (|IID (due to the charging 
effect) is the main modification from the earlier models. Hence, Eq. ( |i~T| ) becomes identical to 
the phase difference equation derived by Bulaevskii et al. (Eq. (11) in Ref. 17) when a = 0. 
Using Eq. ([IT]), we show below that a weak charging effect in BSCCO (i.e., a ~ 0.1 — 0.4) 
can yield significant changes to the phase dynamics. 



III. JOSEPHSON PLASMA DISPERSION RELATION 



We now determine the dispersion relation for the Josephson plasma and the Swihart 
velocity for the collective modes, using linear analysis: <pi,e-i = ff^J-i + ty'ii-x- Here y^,£-i 
describes small fluctuations about </4°iLi describing uniform motion of Josephson vortices in 
the I layer between i th and t — 1 th S layers. </4°i-i is zero in the Meissner state, but in 
general, it depends on a magnetic field,H allowing JPR to be tunned by the field. The effect 
of magnetic field on JPR can be accounted for more accurately via the field dependence of 
J c and via imposing the boundary condition, (cj) /27i)(d(p£ t e-i/dx) = VB, explicitly at both 
x = and junction length) in numerical simulations. 

When the bias current Jb equals the Josephson current in each of the I layers (i.e., 
Jb = J c sin (pf}-i) , we describe the motion of vortices in terms of a uniform motion f\°J_i 
and small perturbation <p'ee-i about (fil-v The uniform phase motion is described by 

AVgL 1 / 1 p dA®\ a 1 ff»gf = 

dx 2 X 2 C V I w\ dt 2 u p dt J X 2 c Vuj 2 dt 2 ' 1 ' 

while small fluctuations (i.e., y^^-i) about y4°-i are described by 
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T Ui ' (°) I I (°) i ' (°) M 

J c [a cos <p\£_ x + s{(p e+1>£ cos (p e+1/ + ¥i-i,e-2 cos ^-1,1-2)] 



8tt 2 dx 2 

+^ jj,e-i + s (Jl+i,e + Jf-\,i-2) 1 (13) 
\u 2 dt 2 ui dt J ~ \ J e+i,e ZJ e,e-i J 1-1^-2) ■ \ L ^) 

Equations (|13D and ([14]) are coupled through the current density Jfi_i, suggesting that these 
equations can be simplified by expressing <^«_i and Jjj>_ x as Fourier series in space for the 
^-direction: ^-1 = ffi 1 ^ 1 and Jj^ = ES'^e'^". k m = mn/(N + l)a 
represents the wavenumber for the collective mode along the z-direction, a = dj + d$, m is 
the mode index, and N represents the number of Josephson junctions in a stack. 

The Josephson plasma mode dispersion relation is determined easily by approximating 
that (fif}_i ~ <fi^ and by combining Eqs. fll3|) and (|l4l) into a single equation as 
d 2 T 
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(k x A c ) .4 m £ m + A m COS 



(0) 



T m w , (15) 



where A m = 1 + 4a sm 2 (k m a/2), B m = {1 + 4[-5/(l + 2S)} sm 2 ( y k m a/2)}~ 1 , and 5 = s/d'. 
Here, we set /3 = for simplicity. From Eq. ([15]) , we obtain the dispersion relation of 
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fc m ) = uj p [(k x X c ) 2 A m B m + A m (cosy? (0) ) t J (16) 

for the collective mode. (• • -) t represents thermal averages. The dispersion relation of Eq. 
(jjl) naturally recovers both purely longitudinal^ and purely transverse plasma excitations^ 
at k x = and at k m — 0, respectively. However, there are notable differences between our 
result of Eq. (|16|) and the results from other models.0il Figure 3 illustrates the difference 
between the dispersion relation of our model and that of the magnetic induction modelll! 
(Fig. 3(a)) and that of the charging effect modeJll (Fig. 3(b)). 

The changes in the dispersion relation due to the charging effect increase the character- 
istic velocity of the collective mode. The group velocity for the electromagnetic waves in 
these LJJ is easily determined from Eq. (|i~6|) by evaluating 



^ = ^ k x \\ A m B m (17) 
dk x uj 

within the linearized model. This group velocity, asymptotically (i.e., as k x — > oo), leads to 
the Swihart velocity 

1 + 4asin (k m a/2) 
.! + 4 (irfs) sin 2 (fc m a/2)J 

the effective maximum velocity for the collective mode m. Here c Q = c/y/e. Equation 

(|l8l) recovers the result of the magnetic induction model (c^f 7 ) when a = (i.e., c^f 7 = 

c m {a = 0)), indicating that the charging effect yields the mode-dependent enhancement 

of the Swihart velocity from c^f 7 . For example, the Swihart velocity is not enhanced for 

the m = 1 mode (i.e., C\ = cf 11 ), but it is enhanced for the m = N mode (i.e., cjy = 

(1 + 4a) 1/2 c^ 7 ). This enhancement reflects the increase in the coupling strength between 

the junctions due to the charging effect and indicates that the threshold velocity v t h (=cjv) for 

emitting Cherenkov radiation^ (i.e., non-Josephson emission) is also increased. For example, 

Vth = c^ 1 when a = 0, but v t h = 1.34cjJ/ 7 when a = 0.2. Evidence, indicating the need 



to account for the charging effect, may be also found in the I-V data for BSCCO. Recent 
analysis of the I-V data in the dense vortex regime (i.e., B ^> H Q ) indicates that a better 
agreement between the predicted and observed position of the kinks can be obtained if the 
Swihart velocity for m > 1 is slightly larger than c^f 7 .§ This suggests that accounting for the 
charging effect is important for quantitative understanding of the kinks in the I-V curves. 

In Fig. 4, we compare the Josephson plasma mode dispersion for (a) Nb-Al/AlO^-Nb 
multilayers and (b) BSCCO in the Meissner state (i.e., {<p^)t ~ 0), using of Eq. fll6[) . To 
illustrate the difference between the dispersion of collective mode for these two systems, 
we use the experimental values for the parameters a (i.e., charging effect strength) and S 
(i.e., induction coupling strength). For the spectrum corresponding to the Nb-Al/A10 x -Nb 
multilayers (Fig. 4(a)), we chose a = 0.0 and S ~ —0.47 (assuming A ~ 900A, dj ~ 20A, 
and ds ~ 3OA).0 Here we chose a = since the charging effect is negligible when d$ is 
much larger than an atomic length. For the spectrum corresponding to BSCCO (Fig. 4(b)), 
we chose a = 0.2 and S ~ -0.49999 (since A ~ 1500A, dj ~ 151, and d s ~ 3A)H Here 
a = 0.2 is chosen. There are two notable differences between Figs. 4(a) and 4(b). First, 
due to a stronger inductive coupling (i.e., S = —0.49999 versus —0.47), the frequency uj/ujp, 
for a fixed k x X c , near k m a = in Fig. 4(b), decreases more sharply with k m a than that in 
Fig. 4(a). Second, due to the charging effect (i.e., a = 0.2 versus 0.0), the collective mode 
frequency for k x X c = shows a dispersion as a function of k m in Fig. 4(b), indicating purely 
longitudinal excitations, while no dispersion is shown in Fig. 4(a), indicating the absence of 
these excitations. Note that the effect of finite, but small, j3 is negligible, here. 

IV. STABILITY OF COLLECTIVE MODES 

In this section, we discuss the stability of uniform motion of collective modes (i.e., moving 
JVL) shown in Fig. 2 against small fluctuations. The structure of the moving JVLs, driven 
by a bias current, evolves as a function of its velocities.0 This evolution can be easily 
understood in terms of the stable-unstable transition for the collective modes. 
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We now carry out the linear analysis and determine the condition for maintaining stable 
uniform motion (i.e., the condition for bound oscillations of ^^-1) by computing the ve- 
locities at which the driven collective modes are stable. Here, instability of uniform motion 
arises when the amplitude of fluctuations grows exponentially as the collective modes propa- 
gate along the junction layers. Similar analysis, not including the charging effect, have been 
carried out to investigate the stability of moving JVL against lattice deformation.0 Also the 
effects of quantum and thermal fluctuations have been studied.i We note that accounting for 
either the dynamic phase transitionil induced by the lattice displacements or the fluctuation 
effects are beyond the scope of the present analysis. 

We proceed the analysis writing the spatial and the temporal dependence of phase fluc- 
tuations (i.e., <Pg£_i) of Eqs. ( |13D and ([[4]) in Fourier space for the z-direction. Combining 
Eqs. ( pT3[) and (|T4l) in Fourier space, we obtain 

\\V d 2 T m A m d 2 T m ( (3 dT n 



d'C m dx 2 uj 2 dt 2 \ujp dt 



+ cos^T m =0 (19) 



where C m = 1 + 2S cos k m a. The uniform motion of the phase locked mode (pj® with the 
wavenumber k m is given by (p$ = u m t + k x x + y? ,m3^ where to m = (2n / 1 <^^)V m A m is the 
Josephson frequency, V m is the average voltage, and ip ,m is a mode dependent constant .0 
Here, the induced field contribution to ipffl from the Josephson effect is neglected. This 
contribution is negligible when the magnetic vortices are in every I layers, as the case for 
B ~ H Q . u m /k x is the velocity of the collective mode m. We transform Eq. ([19]) into 
a familiar Mathieu equation in the following two steps: first, make a change of variables 
from (x,t) to Cm(— )> an d second, let T m = T m e~ VmC - m l 2 . Here T m = (3u) m uj p C m /£l m and 
fl m = AmCmU^ — {k x \ c ) 2 (V I d')uj 2 . The stability condition is determined, solving 

+ [5 T m + C cos (A = , (20) 

where 5 m = — and fj m = C m uj 2 /fl m are the mode dependent (i.e., k m ) parametric 
constants. Solutions of Eq. ( p0|) exhibit instabilityHHS for certain values of 5 m and fj^, 
indicating that the collective mode becomes unstable against small fluctuations. We deter- 
mine the stability condition, finding (5 m , fj^) at which all solutions of Eq. ( p0|) are bounded. 

11 



Note that a similar parametric instability (in the 5^ — r/^ space) occurs both in the magnetic 
induction mo deli (i.e., a = 0) and in the charging effect mo delH (i.e., k x X c = 0). 

For finding the stability condition for the collective modes, it is useful to determine, first, 
the stability diagram of the Mathieu equation 

fL + [5 + V cosQT = 0, (21) 

and then, find the values of (<5^, fj^) corresponding to the stable region of this diagram. The 
boundary curves separating the region of bound (stable) and unbound (unstable) solutions 
can be obtained easily, solving Eq. (|2l]) numerically following the procedure outlined in Ref. 
35. The boundary curves for the periodic oscillations with the period 2ir (i.e., ( = 27r) and 
47r (i.e., ( = Att) are obtained, imposing that the determinant £ n for n = oo, derived from 
Eq. (HD writing T = En=-oo C n e in< for ( = 2vr and T = En=-oo d n e in ^ 2 for ( = 4vr, is zero 
(i.e., £ n = 0).l! Here £ n is the determinant of a (2n + 1) x (2n + 1) matrix for a periodic 
solution with the period 2ir (or a 2n x 2n matrix for a periodic solution with the period Air). 
The determinant £ n can be computed using the recursion relationil 

£n+2 = (1 - 7n+27n+l)£n+l ~ 7n+27n+l(l ~ ln+2ln+l)£n + 7n+27n+l7n^n-l (22) 

where £ = 1, £\ = 1 - 27071, £ 2 = (1 - 7172) 2 - 27071(1 - 7172), and j n = 77/ [2(6 - n 2 )] 
for a solution with the period-27r and £\ = 1 — 7 2 , £ 2 = (1 — 7172) 2 — 7i, ^3 = (1 ~ 7i72 — 
727s) 2 — 7i(l — 727s) 2 , and 7 n = 2rj/[A5 — {2n — l) 2 ] for a solution with the period-47r. 

The stability diagram for Eq. fl2~I|) is shown in Fig. 5. The unstable regions, where 
at least one solution is unbounded, are shaded, and the stable regions, where all solutions 
are bounded, are not shaded. The boundary curves separating these regions are periodic 
solutions with period 2ir (dashed lines) and 4n (solid lines). These curves are obtained by 
calculating £ n = for n = 200. The filled squares represent the values of (5^, 77^) satisfying 
the stability condition. Here we set C = (i.e., (3 = 0) since 5^ oc (3 2 and the terms of this 
order 0{j3 2 ) have been neglected due to small (3. For 5^ = 0, the following values satisfy the 
stability condition: the values shown in Fig. 5 are < fj^ < 0.4540, ff^ ~3.7898, 10.6516, 
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and 20.9637, and the values not shown in Fig. 5 are fj^ ^34.7142, 51.9022, 72.52784, 96.5910, 
124.0918 • • -. For a large fj^ satisfying the stability condition, we assume that Q m « since 
&"m — ► as /^m — * oo. In this case, the velocity for uniform motion is given by 

u m (l + 2S\ l/2 

T^ XcUJp {j~r) ' (23) 

rv x \ ,Ai m i^ m / 

indicating that the presence of the charging effect yields the mode dependent modification to 
the stability condition. For example, u>i/k x ~ \ c u p for m — 1 (i.e., rectangular lattice) but 
u N /k x « A c cj p [(l + 25)/(l-25)(l + 4a)] 1 / 2 for m = N (i.e., triangular lattice). The velocity 
for the out-of-phase modes is reduced from the predicted value of the magnetic induction 
model (i.e., a = 0). Equation ( p3[) indicates that moving Josephson vortices in a periodic 
array evolve from one stable mode to another as the vortex velocity increases. For example, 
as the vortex velocity exceeds Ojsr/k x , but less than oj^_i/k x , the moving triangular lattice 
(m = N) becomes unstable and the m = N — 1 mode becomes stable. 

V. MULTIPLE SUB-BRANCHING OF SUPERCURRENT 

In the resistive state, the supercurrent branch splits into multiple sub-branches as the 
bias current exceeds the Josephson current .i'iil Note that these supercurrent sub-branches 
differ from the observed multiple quasiparticle branchesB0Ei in the I-V data for L JJ stacks. 
This supercurrent sub-branching phenomenon, which appears clearly in the non-linear I-V 
regime, is attributed to the motion of Josephson vortices, but its origin is not understood 
clearly. Microwave induced voltage stepsS and geometric resonancei! are considered as 
other mechanisms, but we do not discuss them here. Instead, we argue that the splitting 
of the supercurrent branch is indeed due to the collective motion of Josephson vortices 
examining the low bias current regime where the I-V characteristics is linear. An analytic 
calculation is more tractable in this regime. Here, we illustrate qualitatively, rather than 
quantitatively, how the charging effect modifies the supercurrent sub-branches since the 
particle-hole imbalance effectQ neglected in this study may also need to be included for a 
quantitative comparison with the I-V data. 
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The current- voltage relations in the resistive state is obtained easily by noting that an AC 
voltage ripple with the Josephson frequency uj£/-i, in addition to the DC voltage, appears 
across the junction when a bias current (Jb), greater than the critical current, is applied. 
This AC voltage ripple is due to the electron-pair tunneling current across the junction. 
Using the modified AC Josephson relation of Eq. (g), the time dependence of the phase 
difference between £ th and £ — 1 th S layers can be written asS 

(fe,e-i{t) ~ (fe,e-i{0) + cde,e-it + tr sina>^_it (24) 

where = (2ir / <j) )[(V^e-i) + a{(V i+ i, e ) - 2{Ve,e-i) + (Ve_i/- 2 ))] is the Josephson fre- 

quency, (V^e-x) is the DC voltage (i.e., time averaged) across the superconductor layers £ 
and £ — 1. V £ ^_ 1 is the amplitude of the AC voltage ripple. This time dependent phase 
difference of Eq. ( f2~4"|) yields a DC critical current responsel! of 

J c siny2 V -i(t) = -Ji (p¥y=±) sin^_i(0) (25) 

\2tt Wt,£-1 J 

across the £ th and £ — 1 th S layers, indicating that the junction becomes resistive when 
Jb exceeds the DC critical current. Here Ji(x) is the first order Bessel function of the first 
kind. Equation fl2"5| ) indicates that the current J^_i = Jb — ^c s i n W-i(0 between two 
adjacent S layers is not uniform along ^-direction, even though a uniform bias current is 
applied. Hence, in this resistive state, we may reduce Eq. (PD to 



1 d 2 A e + £dAA a d 2 E e 1 



u 2 dt 2 Up dt J u) 2 dt 2 J c . 



d'Je,£-l + s(J£+l,£ + Jt-l,£-2j 



(26) 



Here, we neglected the spatial dependence (i.e., x variation) of (ftj-i for simplicity. To 
explicitly express the I-V relation for each collective mode, we now rewrite Eq. (p6|) in 
Fourier space for the z-direction as 

•Am 9 (p m P 9ip m J m /97\ 

u 2 dt 2 u~~df ~ T c ' 1 j 

Note that the first and the second term on the left hand side of Eq. Q27D represent the 
capacitive and the resistive contribution of the junction, respectively. 
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We now average Eq. ( p7[) over time. Since AC Josephson tunneling leads to a small 
voltage oscillation about the DC voltage, a further simplification of Eq. (E7[) can be made. 
The modified AC Josephson relation of dip m /dt = (2ir/(j) )V m Am indicates that the capaci- 
tive contribution vanishes when it is averaged over time (i.e., (d 2 (p m /dt 2 ) oc d(V m )/dt ~ 0). 
This simplification leads to the current-voltage relation of 

< 28 > 

Here (• • •) denotes the time average, and J = (J m ) ~ Jb — J c (sin<^). Since the collective 
modes for m = 1,2, ---N are identical to the modes for m = 2N+1, 2N, •••iV+2, respectively, 
the number of sub-branches is the same as the number of junctions (i.e., N) in the stack. 

In Fig. 6, we plot the I-V relation of Eq. (^) for (a) a = 0.0 and (b) 0.1 to illustrate 
the effect of weak, but non-zero, charging effect. For clarity, we plot the curves for only 
the three collective modes (i.e., m = 1, N/2, and N) corresponding to the modes shown in 
Fig. 2. These I-V curves reveal two interesting points. First, the supercurrent splits into 
iV sub-branches, each corresponding to the collective mode, when the LJJ stack are in the 
resistive state. The m = 1 mode represents the high velocity mode (i.e., rectangular lattice), 
while the m = N mode represents the low velocity mode (i.e., triangular lattice). Second, 
these N sub-branches appear as a single curve when the charging effect is absent (i.e., a = 0, 
see Fig. 6(a)), but they spread out when this effect is present (i.e., a ^ 0, see Fig. 6(b)), 
suggesting that this can be used as another way to measure the charging effect. Since the 
width of this spread is related to the strength of the charging effect (a), identification of 
each sub-branches is feasible at a low bias current. When a is small, as in BSCCO (i.e., 
a ~ 0.1 — 0.4), observing the branch splitting in the linear I-V regime may be difficult but is 
still possible. The main difficulty is in observing the high velocity branches (i.e., m ~ 0(1)). 
To observe these branches, a magnetic field, stronger than B ~ H a , may be needed because 
of their stability conditions. Note that the appearance of these high velocity branches is 
expected when the interaction between the vortices is increased by the field,! suggesting that 
a complete sub-branch structure may be more easily obtained from the I-V characteristics of 
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a LJJ stack with increasing microwave irradiation power (i.e., AC magnetic fields) These 
results are consistent with the dataSHi exhibiting supercurrent branch splitting. 

VI. SUMMARY AND CONCLUSION 

In summary, we investigated the collective phase dynamics in the conventional LJJ stacks 
and in layered superconductors, using a theoretical model which accounts for both the mag- 
netic induction effect and the charging effect .Bill These two coupling mechanisms are equally 
important in the intrinsic LJJ (e.g. BSCCO) due to the atomic length thick S layers. We 
showed that the collective phase dynamics in an intrinsic LJJ stack is modified from those 
in a conventional LJJ stack in two important ways, (i) The dispersion of Josephson plasma 
mode for BSCCO is significantly changed from the Nb-Al/AlO^-Nb multilayers. Conse- 
quently, the Swihart velocity and the velocity of stable uniform motion for the out-of-phase 
collective modes in BSCCO increases and decreases, respectively, from the results of the 
magnetic induction model due to the presence of the charging effect, (ii) The supercurrent 
sub-branching in the resistive state is consistent with collective motion of Josephson vortices. 
The width of spread of these supercurrent sub-branches in the linear I-V regime depends 
on the strength of the charging effect. These results are consistent with the experimental 
data and illustrate the importance of accounting for the charging effect in BSCCO, even 
though its strength is weak (a ~ 0.1 — 0.4). They also suggest that our model is useful for 
understanding the experimental data for JPR, non- Josephson emission, and the I-V char- 
acteristics in the resistive state. Since many applications of intrinsic LJJ stacks as high 
frequency devices exploit collective dynamics of Josephson vortices, these results indicate 
that our model is useful for future technological applications involving intrinsic LJJ stacks. 
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FIGURE CAPTIONS 



Figure 1. A stack of LJJ is shown schematically as alternating superconducting (S) and 
insulating (I) layers with thickness d$ and di, respectively. L x denotes the dimension in 
x— direction. The magnetic field B is applied in the plane of tunnel barriers and the bias 
current Ib is applied along the vertical stack. 

Figure 2. Mutual phase-locking of Josephson vortices (filled ovals) with the wave number 
k m = mir/(N +l)a is schematically illustrated. N is the number of LJJ in the stack and a = 
ds + di. For the in-phase mode (m = 1), the Josephson vortices form a rectangular lattice, 
but for the out-of-phase mode (m — N), they form a triangular lattice. The dotted lines are 
a guide to eyes for phase-locking, and the arrows indicate the direction of propagation. 

Figure 3. The dispersion relation for (a) the longitudinal plasma excitations at k x = 
and (b) the transverse plasma excitations at k m = are plotted to illustrate the difference 
between (a) our model and the magnetic induction model, and (b) our model and the 
charging effect model. Here a = 0.2 and S = —0.49999 are chosen. 

Figure 4. The plasmon dispersion in the Meissner state (i.e., </4°-i = 0) * s plotted as 
functions of k m a and k x \ c for the parameters corresponding to (a) the Nb-Al/AlO^-Nb 
multilayers {a = 0.0, S = -0.47) and (b) BiaSrsCaCuaOg+j, {a = 0.2, S = -0.49999). 

Figure 5. Stability diagram of Mathieu's equation. The unstable regions and the stable 
regions are shaded and not shaded, respectively. The periodic solutions of period 2ir (dashed 
lines) and 47r (solid lines) represent the boundary between the stable and unstable regions. 
The filled squares represent the values of (5 m , fj m ) satisfying the stability condition. 

Figure 6. The I-V curves for the supercurrent branch in the resistive state are plotted for 
(a) a = 0.0 and (b) 0.1 to illustrate the effect of nonzero a (i.e., charging effect). Here, only 
three curves corresponding to the collective modes shown in Fig. 2 are plotted for clarity. 
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